function [ alpha , b, errors, SVMcrits ] = SMO( K,t,C )
%SMO Sequential Minimal optimization
%
% INPUT:
% K			-- Kernel matrix (n x n)
% t			-- Label (double)
% C			-- Regularization constant
%
% OUTPUT:
% alpha     -- alpha coefficients (allow the computation of the bias and
%               core of the classifier)
% b         -- bias
% errors    -- (optional) error at each iteration
%
% USAGE:
% [ alpha , b ] = SMO(K,t,C)
n=length(t);

%% Initialisation
alpha=zeros(n,1);
f=-t;
I_low=find(t<0);
I_up=find(t>0);

nb_iter=0;
last_error=Inf;
errors=[];
SVMcrits=[];
while true
    nb_iter=nb_iter+1;
    
    %% Main loop
    [i, j] = select_pair(f,I_low,I_up);
    if j==-1
        break
    end
    if mod(nb_iter,20)==0
        if f(i)-f(j)==last_error
            fprintf('Not Converging :-/, final error : %f \n',last_error);
            break;
        end
        last_error=f(i)-f(j);
        errors=[errors,last_error];
        %SVM criterion
        %SVMcrit=sum(alpha)-1/2*(alpha.*t)'*K*(alpha.*t);
        %SVMcrits=[SVMcrits,SVMcrit];
    end
    sigma=t(i)*t(j);
    w=alpha(i)+sigma*alpha(j);
    L=max(0,sigma*w-(sigma>0)*C);
    H=min(C,sigma*w+(sigma<0)*C);
    eta=2-2*K(i,j);% diag coeffs replaced by 1
    if eta>1e-15
        alpha_j=alpha(j)+t(j)*(f(i)-f(j))/eta;
        alpha_j=min(max(alpha_j,L),H);
    else
        % Compute Phi_L and Phi_H
        L_i=w-sigma*L;
        Phi_L=0.5*(K(i,i)*L_i^2+K(j,j)*L^2)+sigma*K(i,j)*L_i*L...
            +t(i)*L_i*v_i+t(j)*L*v_j-L_i-L;
        H_i=w-sigma*H;
        Phi_H=0.5*(K(i,i)*H_i^2+K(j,j)*H^2)+sigma*K(i,j)*H_i*H...
            +t(i)*H_i*v_i+t(j)*H*v_j-H_i-H;
        % Update alpha_j
        if Phi_L > Phi_H
            alpha_j=H;
        else
            alpha_j=L;
        end
    end
    alpha_i=alpha(i)+sigma*(alpha(j)-alpha_j);
    f=f+t(i)*(alpha_i-alpha(i))*K(:,i)+t(j)*(alpha_j-alpha(j))*K(:,j);
    alpha(i)=alpha_i;
    alpha(j)=alpha_j;
    %% Update the I_up and I_low sets
    I_up=I_up(I_up~=i & I_up~=j);
    I_low=I_low(I_low~=i & I_low~=j);
    if alpha_i<1e-8
        if t(i)>0
            I_up=[I_up;i];
        else
            I_low=[I_low;i];
        end
    elseif alpha_i>C-1e-8
        if t(i)<0
            I_up=[I_up;i];
        else
            I_low=[I_low;i];
        end
    else
        I_up=[I_up;i];
        I_low=[I_low;i];
    end
    if alpha_j<1e-8
        if t(j)>0
            I_up=[I_up;j];
        else
            I_low=[I_low;j];
        end
    elseif alpha_j>C-1e-8
        if t(j)<0
            I_up=[I_up;j];
        else
            I_low=[I_low;j];
        end
    else
        I_up=[I_up;j];
        I_low=[I_low;j];
    end
end

fprintf('Nb iterations : %i \n',nb_iter);

b=min(f(I_up));


end

function [i_low, i_up] = select_pair(f,I_low,I_up)
theta=1e-8;
%% Compute i_low i_up
[~,i_up]=min(f(I_up));
i_up=I_up(i_up);
[~,i_low]=max(f(I_low));
i_low=I_low(i_low);
%% Check for optimality
if f(i_low) <= f(i_up)+2*theta;
    i_low=-1;
    i_up=-1;
end
end
